//888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//	File		 :	ParseMSF.cpp
//	Contents	 :	Source file for class CSequence
//	Creation Date:	22/12/1999
//	Last Modified:	23/12/1999	lg
//	tabs width	 :	every 4 characters
//
//888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#include <stdexcept>
#include <algorithm>
#include <vector>
#include <string>
#include <cassert>
#include "sequence.h"
#include "std_helper.h"



// private functions manipulating sequences
namespace
{
unsigned GetNamesMaxLen(const CSequence& seq);
void RightJustifyNames(CSequence& seq);
void LeftJustifyNames(CSequence& seq);
void CentreJustifyNames(CSequence& seq);
void ConstructLookupTable(char* p, bool allowLower, const std::string& invalid);
const char isvalid(const char* LookupTable, char c);
//_______________________________________________________________________________________
//
//	ElideSequence
//
//		removes a section in the supplied sequence
//			from begin until end (exclusive aka STL)
//		cntResidues is an array containing the number of residues to be elided
//_______________________________________________________________________________________
void ElideSequence(CSequence& seq, unsigned begin, unsigned end, unsigned* cntResidues, char* lookupTable);
}
#include "sequence.cxx"


//_______________________________________________________________________________________
//
//	GetNamesMaxLen
//
//		Get the longest length of names (including extras)
//_______________________________________________________________________________________
unsigned CSequence::GetNamesMaxLen() const
{
	unsigned maxLenNames = 0;
	for (int i = Count() - 1; i >= 0; --i)
		maxLenNames = std::max(	static_cast<unsigned>(Names[i].length()), 
					maxLenNames);
	for (int i = ExtraCount() - 1; i >= 0; --i)
		maxLenNames = std::max(static_cast<unsigned>(ExtraNames[i].length()), 
					maxLenNames);

	// not really necessary because should cause CreateSequenceName first
	// consensus 80%:
	// 01234567890123
	return(std::max(14U, maxLenNames));
}


//_______________________________________________________________________________________
//
//	Left/Right/Cener JustifyNames
//
//		justify sequence nameswith whitespace
//		Last modified: lg 23/12/1999
//
//_______________________________________________________________________________________
void CSequence::JustifyNames(COptions::JUSTIFY_FLAG f)
{
	switch (f)
	{
		case COptions::JUSTIFY_L:
			LeftJustifyNames(*this);
			break;
		case COptions::JUSTIFY_C:
			CentreJustifyNames(*this);
			break;
		default:
		//case COptions::JUSTIFY_R:
			RightJustifyNames(*this);
			break;
	}
	
}	


//_______________________________________________________________________________________
//
//	CheckLengthsEq
//
//		returns false if any of the sequences have difference data lengths from the first
//		Last modified: lg 07/04/2000
//
//_______________________________________________________________________________________
bool CSequence::CheckLengthsEq() const
{
	unsigned firstSeqLen = Data[0].length();
	for (int i = Data.size() - 1; i > 0; --i)
		if (Data[i].length() != firstSeqLen)
			return false;
	for (int i = ExtraData.size() - 1; i > 0; --i)
		if (ExtraData[i].length() != firstSeqLen)
			return false;
	return true;
}

//_______________________________________________________________________________________
//
//	CheckDuplicateNamesCnt
//
//		returns the number of sequences with duplicate names
//		Last modified: lg 07/04/2000
//
//_______________________________________________________________________________________
unsigned CSequence::CheckDuplicateNamesCnt() const
{
	if (!Count())
		return 0;
	std::vector<std::string> tempNames(Names);
	std::sort(tempNames.begin(), tempNames.end());
	std::vector<std::string>::iterator pDuplicateNames
						= std::unique(tempNames.begin(), tempNames.end());
	return tempNames.end() - pDuplicateNames;
}

//_______________________________________________________________________________________
//
//	CheckEmptySequences()
//		returns first empty index
//	RemoveEmptySequences()
//
//		Checks for and removes empty sequences
//		Last modified: lg 07/04/2000
//
//_______________________________________________________________________________________
int CSequence::CheckEmptyData() const
{
	for (int i = Data.size() - 1; i >= 0; --i)
		if (!Data[i].length() == 0)
			return i;
	return -1;
}
void CSequence::RemoveEmptyData()
{
	// remove sequences with no data backwards so that indices remain valid
	for (int i = Data.size() - 1; i >= 0; --i)
		if (0 == Data[i].length())
		{
			Names.erase(Names.begin() + i);
			Data.erase(Data.begin() + i);
		}

}


//_______________________________________________________________________________________
//
//	Remove
//
//		removes data and name at an index if valid
//		(stretches with > thresholdEmpty # of gaps
//			& of > minimumLen length)
//		Last modified: lg 06/4/2001
//
//_______________________________________________________________________________________
void	CSequence::Remove(unsigned index)
{
	assert(IS_VALID());
	if (index >= Count())
		return;

	Names.erase(Names.begin() + index);
	Data.erase(Data.begin() + index);
}

//_______________________________________________________________________________________
//
//	Excise
//
//		Cuts out a stretch of residues from s to e not inclusively, i.e. (s, e]
//		Last modified: lg 06/4/2001
//
//_______________________________________________________________________________________
void CSequence::Excise(unsigned s, unsigned e)
{
	assert(s < Len && e <= Len && s < e);
	for (unsigned i = 0; i < Data.size(); ++i)
		Data[i].erase(s, e - s);
	for (unsigned i = 0; i < ExtraData.size(); ++i)
		ExtraData[i].erase(s, e - s);
	Len = Data[0].length();
}

//_______________________________________________________________________________________
//
//	RemoveGaps
//
//		Cuts out a regions where all the sequences in a column consist of gaps
//		Last modified: lg 06/4/2001
//
//_______________________________________________________________________________________
void CSequence::RemoveGaps(const char* loopupTable)
{
	// work from end to beginning to preserve indices as
	// the sequences are cut shorter
	for (int endPos = Len - 1; endPos >= 0; --endPos)
	{
		// find the end of a gapped region
		// check if each position is all gaps
		unsigned cntResidue = 0;
		for (int j = Count()-1; j >= 0 && !cntResidue; --j)
			cntResidue += isvalid(loopupTable, Data[j][endPos]);
		
		if (!cntResidue)
		{
			// now find the start of this gapped region
			// actually finds start - 1
			int i;
			for (i = endPos - 1; i >= 0; --i)
			{
				unsigned cntResidue = 0;
				for (int j = Count()-1; j >= 0 && !cntResidue; --j)
					cntResidue += isvalid(loopupTable, Data[j][i]);
				
				if (cntResidue)
					break;
			}
			
			Excise(i+1, endPos+1);
			
			// start seaching again before this sparse region
			endPos = i-1;
		}
	}
	
	assert(Len == Data[0].length());
	assert(IS_VALID());
}

//_______________________________________________________________________________________
//
//	ElideSparseRegions
//
//		removes sparse regions
//		(stretches with > thresholdEmpty # of gaps
//			& of > minimumLen length)
//		Last modified: lg 06/4/2001
//
//_______________________________________________________________________________________
void CSequence::ElideSparseRegions(unsigned thresholdGaps, unsigned minimumLen,
								bool allowLowerCase, const std::string& invalidStr)
{
	char lT[256];
	ConstructLookupTable(lT, allowLowerCase, invalidStr);

	RemoveGaps(lT);
		
	assert(minimumLen >= 4);
	unsigned thresholdResidue = Count() - thresholdGaps;
	unsigned* cntResidues = new unsigned[Count()];


	
	// work from end to beginning to preserve indices as
	// the sequences are cut shorter
	for (int endPos = Len - 1; endPos >= 0; --endPos)
	{
		// find the end of a sparse region
		// check if each position is mostly blank
 		unsigned cntResidue = 0;
		for (int j = Count()-1; j >= 0 ; --j)
		{
			cntResidue +=
			cntResidues[j] = isvalid(lT, Data[j][endPos]);
		}

		// now find the start of this sparse region
		// actually finds start - 1
		// if residues < threshold, still in sparse
		if (cntResidue <= thresholdResidue)
		{
			int i;
			for (i = endPos - 1; i >= 0; --i)
			{
				unsigned cntResidue = 0;
				for (int j = Count() - 1; j >= 0 ; --j)
				{
					if (isvalid(lT, Data[j][i]))
					{
						++cntResidues[j];
						++cntResidue;
					}
				}
				if (cntResidue > thresholdResidue)
					break;
			}

			// elide if sparse region > minimum length
			if (endPos - i >= (int)minimumLen)
			{
				// end pos points to 1 past sparse region, wind back
				if (i > 0)
				{
					for (int j = Count() - 1; j >= 0 ; --j)
						if (isvalid(lT, Data[j][i]))
							--cntResidues[j];
				}
				ElideSequence(*this, i+1, endPos+1, cntResidues, lT);
			}

			// start seaching again before this sparse region
			endPos = i-1;
		}
	}

	// reset Len because sequences have been elided shorter
	Len = Data[0].length();
	assert(IS_VALID());
	delete[] cntResidues;
}



//_______________________________________________________________________________________
//
//	CreateConsensusName
//
//		Sets the name of the first ExtraData line to "Consensus"
//		Last modified: lg 11/04/2000
//
//_______________________________________________________________________________________
void CSequence::CreateConsensusName(bool useThreshold, unsigned threshold)
{
	assert (ExtraCount() >= 1);
	if (ExtraCount() < 1)
		return;

	ExtraNames[0] = "Consensus";
	if (useThreshold)
	{
		char buf[20];
		sprintf(buf, "/%d%%", threshold / 10);
		ExtraNames[0] += buf;
	}
}


//_______________________________________________________________________________________
//
//	ChkValidity
//
//		returns a vector of errant characters within each sequence in the alignment
//		Last modified: lg 11/04/2000
//
//_______________________________________________________________________________________
bool CSequence::ChkValidity(std::vector<CSequenceInvalid>& invalidSequences,
							bool allowLowerCase, bool allowStrangeLinkers, const std::string& invalidResidues) const
{
	for (unsigned j = 0; j < Count(); ++j)
	{
		short cnts[256];
		memset(cnts, 0, sizeof cnts);
		std::string::const_iterator i = Data[j].begin();
		std::string::const_iterator e = Data[j].end();

		while (i != e)
			++cnts[(unsigned char)*i++];
		CSequenceInvalid invalid;

		// check for whitespace
		unsigned whitespace =	cnts['\n'] +
							cnts['\r'] +
							cnts[' ']  +
							cnts['\t'];
		invalid.Whitespace = whitespace != 0;

		// check for non-printing
		unsigned nonprinting = 0;
		for (int i = 1; i <= 32; ++i)
			nonprinting += cnts[i];
		nonprinting -= whitespace;
		invalid.Nonprinting = nonprinting > 0;


		// check for invalid residues
		std::string invalidRes = invalidResidues;
		for (unsigned i = 0; i < invalidRes.length(); ++i)
			if (cnts[invalidRes[i]])
				invalid.InvalidRes += invalidRes[i];

		// check for lower case
		if (!allowLowerCase)
		{
			for (char i = 'a'; i <= 'z'; ++i)
				if (cnts[i])
					invalid.LowercaseRes += i;
		}
		else
		// check for lower case invalid residues
		{
			for (unsigned i = 0; i < invalidRes.length(); ++i)
				if (cnts[invalidRes[i] - 'A' + 'a'])
					invalid.InvalidRes += invalidRes[i];
		}

		// check for strange linkers
		if (!allowStrangeLinkers)
		{
			for (unsigned i = 33; i < 'A'; ++i)
				if(cnts[i])
					invalid.NonStdLinkers += (char)i;
			for (unsigned i = 'Z' + 1; i < 'a'; ++i)
				if(cnts[i])
					invalid.NonStdLinkers += (char)i;
			for (unsigned i = 'z' + 1; i <= 255; ++i)
				if(cnts[i])
					invalid.NonStdLinkers += (char)i;


			std::string std_Linkers("~-.");
			RemoveCharsFromString(invalid.NonStdLinkers, std_Linkers);
		}

		if (invalid.isInvalid())
		{
			invalid.Name = Names[j];
			invalidSequences.push_back(invalid);
		}
	}
	return invalidSequences.size() == 0;
}


//_______________________________________________________________________________________
//
//	ToUpper_and_Substitute
//
//		forces letters to upper case
//		converts all the characters in "from" to the "to" character
//		perform two sorts of substitution in one pass for speed
//		Last modified: lg 11/04/2000
//
//_______________________________________________________________________________________
void CSequence::ToUpper_and_Substitute(bool toupper, bool substitute, char to, const std::string& from)
{
	// if no changes return
	if (!toupper && !substitute)
		return;
		
	char redirect[256];
//		for (int i = 0; i < 256; ++i)
//			redirect[i] = i;
//		fun optimization!
	unsigned* u = reinterpret_cast<unsigned*>(redirect);
	for (unsigned i = 0, val = 0x03020100; i < 64; val += 0x04040404, ++i)
		u[i] = val;
	if (toupper)
		for (char i = 'a'; i <= 'z'; ++i )
			redirect[i] = (char)(i + 'A' - 'a');
	if (substitute)
	{
		for (int i = from.length() - 1; i >= 0; --i )
			redirect[from[i]] = to;
	}
	for (int i = Count() - 1; i >= 0; --i)
	{
		std::string::iterator p = Data[i].begin();
		std::string::iterator e = Data[i].end();
		while (p != e)
		{
			*p = redirect[(unsigned char)*p];
			++p;
		}
	}
}	


#ifndef NDEBUG
bool CSequence::IS_VALID() const
{
	if (Len ==0 || Len > 32000)
		return false;
	if (Count() > 10000)
		return false;
	if (ExtraData.size() != ExtraCount())
		return false;
	if (Data.size() != Count())
		return false;
	for (unsigned i = 0; i < Count(); ++i)
		if (Data[i].length() != Len)
			return false;
	for (unsigned j = 0; j < ExtraCount(); ++j)
		if (ExtraData[j].length() != Len)
			return false;
	return true;
}
#endif
